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ABSTRACT 


The  differential  equations  for  the  evolution  of  the  fluid  quantities  in  steady 
detonation  waves  are  integrated  numerically  using  a  fourth  order 
Runge-Kutta-Fehlberg  method.  Strong  and  Weak  detonation  waves  are  found  as 
well  as  detonation  failure.  In  the  weak  regime  there  is  a  continuum  of 
propagation  speeds  for  a  given  expansion  loss. 
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INTRODUCTION 


To  avoid  repetition  the  initial  proposal  for  this  grant  is  attached  as  an 
Appendix  to  this  Final  Report.  The  work  in  the  grant  was  to  study  the 
existence  of  multiple  propagation  speeds  and  detonation  failure  for  one 
dimensional  steady  detonation  waves  with  expansion  losses.  The  theoretical 
treatment  was  based  on  the  model  developed  in  Sec. SO  of  the  book  by  Fickett 
and  Davis^.  Only  the  steady  reaction  region  between  the  initiating  shock  front 

and  the  sonic  locus  was  considered.  The  main  fluid  quantities  are  the  Mach 

2  2 

number  squared,  m  =  u  /c  ,  and  the  reaction  progress  variable,  X.  These 
quantities  obey  the  equations 


2(l-m)dm/dx 
md+ym)  (2+(Y~l)m) 


d(^>  A^q-l/dx  2e 

1  ^iqi+^l+ul^  1+Ym 


dx  u 


(1) 


(2) 


and 


2 

c 


2(y-1) 

2  +(Y~l)m 


xi*li  +  uj  + 


The  losses  are  taken  into  account  by  the  parameter  e  as  defined  in  the 
proposal.  In  reality  the  expansion  loss  is  not  a  given  quantity.  It  is 
determined  by  the  complete  three  dimensional  flow  of  the  explosive  charge  and 
the  confining  medium.  This  entails  solving  a  free  boundary  problem  of  great 
complexity.  The  chemical  reaction  rate  R.  is  given  in  Eq.(3)  of  the  proposal. 
It  involves  the  parameters  B  and  G  and  RV  that  appear  in  the  figures  below.  G 
is  merely  the  ratio  of  the  specific  heats,  y,  while  B  and  RV  characterize 
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the  strength  of  the  explosive.  B  was  defined  in  the  proposal,  and  RV  is  given 
2 

by  RV  =  c0/2q.  The  larger  the  heat  of  reaction,  q,  that  is  released,  the 
smaller  B  and  RV  are.  However,  B  increases  with  the  activation  temperature  of 
the  reaction. 


We  are 

in  the  frame 

of  the  reaction  front 

.  The  unreacted  explosive 

fluid 

approaches 

f  rom  x  <  0 

with  a  supersonic 

detonation  velocity 

D 

to  be 

de  termined. 

D  will  be 

obtained  in  terms 

of  the  Mach  number 

squared 

_2.  2 
m0  =  13  !cq 

>1  in  front 

of  the  shock  at 

x  =  0  .  The  flow 

at 

x=0+. 

immediately 

behind  the 

shock,  is  subsonic 

.  The  Mach  number 

squared 

“l  =  ^ul^cl 

)2  at  i  =  0  + 

satisfies  m^  <.  1  and 

is  related  to  mo  by 

the 

usual 

2 

Hugoniot  relations 


2  + 

1  2ymo  ~  Cr-l) 

and 


c 


2 

1 


2  (2ym0-(y-l))(2+(y-l)«0) 

°  m0(y+l)2 


Several  authors  have  considered  this  model  with  different  reaction  rates 

3 

R^  and  expansion  loss  terns  e.  Wecken  considered  the  case  e  =  constant 

and  obtained  partial  results  on  multiple  solutions  and  detonation  failure. 

4  5 

Tsuge  et  al.,  and  Fuj iwara  and  Tsuge  have  integrated  the  nozzle  equations 
numerically  for  a  hydrogen/oxygen  detonation  with  a  good  model  for  the  lateral 
expansion.  They  have  shown  that  sufficiently  great  expansion  losses  lead  to 
detonation  failure.  They  also  obtain  two  detonation  speeds  that  join  at  the 
failure  diameter.  For  small  losses  the  fast  root  goes  to  the  Chapman-Jouguet 
value.  They  conjectured  that  the  slow  root  went  to  Mach  one  as  the  losses 


2 


decreased  to  zero  (charge  diameter — )  “) .  They  were  not  able  to  verify  it  due 
to  the  large  mounts  of  computer  time  required. 

METHOD  OF  SOLUTION 

The  phase  plane  of  Eqs.(l)  and  (2)  is  extraordinarily  complex  even  with 
e  =  constant  as  we  take  it  here.  Let  F(m,  X  ;  m0,  c)  be  the  right  hand  side 
of  Eq.(l)  so  the  curve  F  =  0  (with  mQ  and  e  fixed)  is  the  same  as  the  locus 
of  dm/dX  =0  in  the  m  -  X  plane.  The  singular  points  of  Eqs.(l)  and  (2) 
are  given  by  the  intersections  of  the  line  m  =  1  with  the  curve  F  =  0.  The 
phase  trajectory  of  the  physical  solution  begins  at  the  point  immediately 
behind  the  shock  where  x  «=  0+.  There  X  =  0  while  m  has  the  so  far  unknown 
value  m  =  mj .  The  crucial  point  is  that  the  physical  trajectory  has  to  be 
transonic,  that  is  it  must  cross  the  line  m  =  1  with  a  finite  slope  dm/dX  so 
it  can  be  joined  to  the  rear  flow  behind  the  sonic  point.  Therefore  the 

transonic  trajectory  has  to  cross  the  line  m  =  1  at  a  critical  point.  In 

normal  phase  plane  problems  the  orbits  may  only  go  to  the  critical  points 
asymptotically,  reaching  them  only  after  an  infinite  distance  (or  time).  Here 
the  transonic  trajectories  do  pass  through  the  critical  points  due  to  the 
singularity  at  m  =  1.  The  condition  that  the  solution  must  be  transonic  is 

what  determines  the  value  of  m^  (as  well  as  mQ  which  is  related  to  m^  by  the 

Hugoniot  relation).  The  curve  F  =  0  is  not  very  complicated  qualitatively. 
For  different  values  in  the  physical  range  (y-l)/2y  <  Bj  <  1,  the  curve 
F  =  0  intersects  the  line  m  =  1  a  different  number  of  times  depending  on  the 
values  of  e  and  m0.  In  the  cases  considered  here  they  can  intersect  twice, 
once,  or  not  at  all.  The  difficulty  is  that  the  curve  dm/dX  *  0,  and 

consequently  the  nature  and  position  of  the  critical  points,  as  well  as  the 
phase  trajectories,  depend  on  the  unknown  ®0(or  m^)  that  enters  not  only  as  an 
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initial  value,  but  also  as  a  parameter  in  the  equations. 

A  typical  case  for  a  strong  detonation  solution  is  shown  in  Fig.l  where 
e,  as  well  as  the  parameters  B,  G,  and  RV  are  given  in  the  table  next  to  the 
figure.  The  curves  F  =  0  and  m  =  1  are  shown  in  Fig.l.  The  curve  F  *  0  is 
computed  numerically  using  a  root  finding  algorithm  (Derpar).In  Fig.(l),  the 
curve  F  =  0  is  plotted  for  the  value  mQ  =  39.16  .  The  differential  equations 
for  the  trajectories  are  integrated  using  a  fourth  order  Runge-Kutta-Fehlberg 
method  whose  variable  step  size  allows  important  savings  in  computer  time. 
For  mQ  =  37.16  the  phase  trajectory  starts  at  the  point  (A  =  0,  m«  *  0.3454) 
and  reaches  the  line  m  =  1  with  a  vertical  slope.  This  is  a  trajectory  of  type 
2.  a  choking  trajectory.  For  mQ  =  41.16  the  trajectory  starts  at 

(A  =  0,  mj  =  0.3442),  crosses  the  line  F  =  0  and  acquires  a  negative  slope. 
This  a  trajectory  of  type  1,  a  subsonic  trajectory.  The  transonic  trajectory 
will  have  a  value  of  mo  bracketed  37.16(type  2)  <  (transonic)  <  41.16(type  1) 
(me  *  39.16  is  pretty  close).  It  passes  through  the  critical  point  at  m  =  1, 
F  =  0.  This  is  the  sonic  point.  The  problem  of  finding  mQ  is  not  an  initial 
value  problem,  but  rather  a  sort  of  boundary  value  problem.  Its  form  is 

suitable  for  a  Nonlinear  Shooting  Method**.  We  start  with  a  value  of  mc  and 

get  its  trajectory.  Then  we  change  the  value  of  mo  until  the  type  of  the 

trajectory  changes.  We  have  then  bracketed  the  desired  value  of  mQ.  This 
procedure  is  done  for  each  value  of  c.  In  the  end  we  produce  a  plot  of 
m0  vs.  e  that  shows  how  the  steady  detonation  speed  varies  with  the  loss 
parameter. 

STRONG  SOLUTIONS  THROUGH  HYPERBOLIC  (SADDLE)  CRITICAL  POINTS 

When  the  sonic  point  is  a  hyperbolic  critical  point  (saddle  point)  the 

i 

shooting  method  works  well  because  there  is  a  separatrix  that  separates  the 
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trajectories  of  types  1  and  2  unambiguously.  In  fact,  the  transonic  trajectory 
lies  on  the  siparatrix.  In  all  the  cases  considered  here  the  singular  point  is 
hyperbolic  if  the  curve  F  =  0  has  a  negative  slope  at  the  critical  point 

(whenever  we  speak  of  a  critical  point  it  must  be  remembered  that  its  position 
depends  on  mQ,  so  that  as  we  vary  mQ  to  produce  different  trajectories  we  are 
also  moving  the  critical  point).  The  phase  trajectories  of  strong  detonation 
waves  always  end  at  a  saddle  point.  Fig.l  is  a  typical  example  of  this 
situation.  In  general  there  is  a  range  of  values  of  e  where  for  each  given 
e  one  can  find  a  single  value  of  mo  (and  therefore  a  single  steady 
detonation  speed)  that  yields  a  transonic  trajectory  and  solves  the  problem. 
In  the  range  of  strong  detonations  mfl  is  a  decreasing  function  of  e. 

STRONG  AND  WEAK  SOLUTIONS  AND  DETONATION  FAILURE 

There  is  another  range  of  e  where  there  are  two  detonations  speeds. 

That  is,  for  a  given  e  there  are  two  values  of  mo  that  yield  transonic 

trajectories  through  their  respective  saddle  points.  The  solution  with  the 

larger  value  of  m  =  ms  is  a  strong  detonation  wave  and  the  solution  with  the 

smaller  aw  is  a  weak  detonation  wave.  As  e  increases,  m$  decreases  while  mw 

increases.  An  example  is  given  in  Figs. (2)  and  (3).  The  values  of  e  «  0.06, 

B  =  6.83,  G  =3.0,  and  RV  =  0.102  are  the  same  in  both  figures.  In  the  figures, 

the  notation  F(m  )  =  0  means  that  the  curve  is  shown  for  the  given  value  of  in 
o  o 

with  the  other  parameters  as  given  in  the  accompanying  table.  The  value  of  mQ 
for  the  strong  transonic  trajectory  in  Fig. (2)  is  bracketed  21.57  (type 
2)  <  m^  <  23.57  (type  1).  The  value  ms  ~  22.57  is  bracketed  numerically  much 
more  closely  than  shown  above.  Those  bracketing  values  are  chosen  for  clarity 


in  the  figure. 

The 

strong  detonation 

of  Fig. 

(2) 

has  the  same 

type  of 

bracketing  as 

the 

case  of  Fig. ( 1) , 

that  is 

the 

transonic  value 

o 

S 

o 

m 

5 
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bounded  on  the  low  side  by  a  type  2  solution  and  on  the  high  side  by  a  type  1 
solution.  Fig. (3)  shows  a  representative  weak  detonation  case.  The  transonic 
value  is  bracketed  by  13.0  (type  1)  <  mw  <  15.0  (type  2).  This  bracketing  for 
the  weak  solution  mw  ~  14.06  is  the  reverse  of  the  strong  one  because  the 
lower  m0  value  is  type  1  while  for  the  strong  case  the  lower  value  was  type  2. 
This  is  a  distinguishing  feature  of  the  strong  and  weak  solutions.  Note  that 
the  curve  F  =  0  in  both  Figs. (2)  and  (3)  looks  pretty  much  like  it  did  in 


Fig. 

(i). 

In 

each  case  there  is 

just  one 

saddle 

point.  As  e 

increases. 

the 

weak 

and 

the 

strong  solutions 

approach 

each 

other  until 

they  meet 

a  t 

“o  = 

mdf 

=  17. 

08  when  e  =  j.  «  0 . 

0657  . 

Af  ter 

e  exceeds  this 

value  there 

i  s 

no 

transonic 

trajectory  at  all 

.  This 

point 

is  usually  called  detonation 

failure.  Fig. (4)  shows  what  the  phase  trajectories  and  the  F  =  0  curves  look 
like  just  past  that  point,  say  for  e  -  .08  .  Past  the  detonation  failure 
value  the  trajectories  remain  of  type  1  as  mQ  is  varied,  no  bracketing  occurs, 
and  there  is  no  steady  detonation  solution.  Note  that  there  is  still  a  saddle 
point  in  the  range  of  mQ  shown  but  it  is  not  reached.  One  might  have  thought 
that  detonation  failure  would  only  occur  when  the  F  =  0  curve  stays  so  low 
that  it  never  crosses  the  n  =  1  line  and  there  is  no  critical  point.  The 
actual  limiting  point  of  detonation  failure  does  not  occur  when  the  F  =  0 
curve  is  tangent  to  the  m  =  1  line.  It  is  a  little  more  subtle  than  that. 

The  bracketings  mentioned  above  are  indicative  of  the  behavior  of  the 
phase  trajectories  in  the  regions  of  strong  and  weak  solutions.  In  the  strong 
region  the  type  1  trajectories  move  higher  monotonical ly  as  is  decreased. 
Then  as  one  decreases  mQ  the  trajectory  rises  until  it  becomes  type  2  and  the 
transonic  value  is  bracketed.  In  the  weak  region  the  reverse  happens.  The 
type  1  trajectory  moves  higher  monotonical ly  as  mQ  is  increased  (rather  than 
decreased)  until  it  becomes  type  2  and  brackets  the  transonic  trajectory. 
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Fig. (5)  is  a  schematic  that  helps  to  visualize  what  is  happening.  Detonation 
failure  occurs  when  these  two  behaviors  interfere  with  each  other  so  the 
trajectories  do  not  behave  monotonical ly  with  mQ.  As  we  decrease  m0,  raising 
the  type  1  trajectory  and  expecting  to  reach  a  strong  solution,  the  behavior 
reverses  and  the  type  1  trajectory  begins  to  move  lower  with  further  decreases 
in  mo  as  it  does  in  the  weak  region.  Consequently,  the  critical  point  is  not 
reached  and  no  transonic  trajectory  is  found.  This  is  evident  in  Fig. (4). 
Toward  the  left  the  trajectories  are  ordered  in  the  strong  way  with  45  below 
30  below  15.  Far  to  the  right  they  are  ordered  in  the  weak  way  with  15  below 
30  below  45.  In  the  region  where  the  trajectories  are  crossing  their  order  is 
not  monotonic.  This  is  what  signals  the  proximity  of  the  limiting  point  of 
detonation  failure.  Note  also  how  the  F  =  0  curves  reverse  their  ordering  in 
this  neighborhood. 

The  above  discussion  can  be  stated  in  a  mathematically  more  precise  way. 
In  a  phase  trajectory  of  type  1,  the  function  m(X)  reaches  its  highest  value 
at  the  point  where  the  trajectory  crosses  the  curve  F  =  0  (for  type  2 
trajectories  it  is  more  useful  to  consider  the  function  A(m)).  Call  this 
value  mjj(m0,e)  noting  that  it  is  a  function  of  the  values  of  mQ  and  e  for 
that  trajectory  (unfortunately,  this  function  is  only  known  numerically).  In 
the  strong  region  dm^/dn^  <  0  ,  while  in  the  weak  region  dm^/dn^  >  0  . 
Fig. (6)  shows  schematically  the  types  of  behavior  of  m^  vs.  mQ  for  several 
fixed  values  of  e.  In  Fig. (6a)  e  is  in  the  double  root  region.  Here  m$ 
is  the  strong  root  and  mw  is  the  weak  root.  In  Fig. (6b)  we  are  at  the  point  of 
detonation  failure  where  e  For  that  value  a^  =  ms  and  the 
nh^"o,edf)  curve  is  tangent  to  the  m  =  1  line.  In  Fig. (6c)  we  are  in  the 
failed  region,  e  )  e  j  j.  and  there  is  no  root. 
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CONTINUUM  OF  SOLUTIONS  THROUGH  A  STABLE  NODE 

Another  type  of  solution  can  occur  when  there  are  two  sonic  points.  The 
curve  F  =  0  depends  on  both  e  and  mQ.  As  e  increases  (or  mQ  decreases)  the 
left  side  of  the  F  =  0  curve  bends  down  and  may  go  below  the  m  =  1  line.  Then 
there  are  two  singular  points.  The  one  with  the  larger  value  of  X  is  a 
saddle  but  the  one  with  the  smaller  value  of  X  is  a  stable  node  (sink)  in 
all  the  cases  considered  here.  The  stable  node  does  not  have  a  separatriz 
like  the  saddle  point  does.  Rather  it  has  limiting  curves  that  act  '  .iges 
for  a  continuum  of  type  3  trajectories.  The  type  3  trajectories  are  t;  ionic. 
They  pass  through  the  stable  node  with  a  finite  slope.  The  behavior  ,  'hat 
for  a  given  e  there  is  a  continuous  range  of  m0  values  that  yield  type  3 
trajectories.  Values  of  mo  above  the  upper  limit  of  the  type  3  continuum 
yield  trajectories  of  type  2.  Values  of  mo  below  the  lower  limit  of  the  type  3 
continuum  may  yield  trajectories  of  type  1  or  type  2  depending  on  what  is 
happening  with  the  coexisting  saddle  point.  Therefore,  we  have  two  cases.  The 
case  t2  <  t3  <  t2,  where  the  type  3  trajectories  are  bounded  above  and  below 
by  type  2  trajectories,  and  the  case  tl  <  t3  <t2  where  they  are  bounded  below 
by  type  1  trajectories  and  above  by  type  2.  Fig. (7)  shows  an  example  of  the 
behavior  of  the  coexisting  saddle  point  that  will  lead  to  the  case 
t2  <  t3  <  t2.  The  parameters  B,  G,  and  RV  are  the  same  as  in  Figs.(2)-(4)  and 
e  =  0.0001.  Fig. (7)  shows  the  F  =  0  curve  plotted  for  mo  =  1.035.  The  sink 
is  to  the  left  and  the  coexisting  saddle  point  is  far  to  the  right.  The 
transonic  trajectory  through  the  saddle  point  is  bracketed 
1.03  (type  1)  <  <  1.04  (type  2).  We  call  the  transonic  value  here  ®wj,c 
because  this  is  a  weak  root  (note  that  type  1  brackets  below)  that  lies  below 
the  continuum.  As  we  increase  mQ  the  type  1  trajectory  rises  and  becomes  type 
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2  around  bq  ~  1.04.  We  continue  to  increase  mo  and  the  type  2  trajectory 
keeps  rising,  moving  mostly  to  the  left,  and  we  enter  the  picture  shown  in 
Fig. (8).  Here  the  lowest  type  2  curve  has  mQ  =  1.25.  Increasing  mQ  produces 
higher  type  2  curves,  we  just  show  the  one  for  mQ  =  1.4  .  Increasing  mQ 
further  gets  us  into  the  type  3  curves.  They  are  shown  here  for  mQ  =  1.7,  2.2, 
2.7,  and  3.2  but  there  is  a  continuum  of  these  transonic  solutions.  They  have 
a  characteristic  shape.  They  reach  the  m  =  1  line  with  a  finite  slope  because 
they  are  all  going  through  their  respective  sinks  (recall  that  the  singular 
point  moves  as  we  change  mo) .  Finally  we  pass  the  continuum  of  type  3  and  run 


into  type 

2  curves 

again  near  mo  ~  3.7  . 

This  case  of  1 2  <  1 3  <  t2 

seems 

to 

be  rare 

because 

it  occurs  here  only 

in  a  narrow  range  of  c  and 

doe  s 

not 

occur  at  all  in  some  of  the  other  examples  shown  below  where  we  take  different 
values  for  B  and  RV.  Fig. (9)  shows  the  type  3  continuum  in  the  case 
tl  <  t3  <  t2.  The  trajectory  rises  as  we  increase  mQ.  We  start  with  type  1  at 
mQ  =  2.0.  The  trajectory  rises,  but  is  still  type  1  at  2.5.  Near  mQ  as  2,77 
(not  shown)  the  type  3  trajectories  start.  We  show  type  3  trajectories  for 
mc  =  3.0,  3.5,  4.0,  4.5,  5.0,  and  5.5  but  there  is  a  continuum  of  transonic 
trajectories  here  going  through  their  respective  singular  points.  The 
continuum  ends  near  mQ  ~  5.5  and  we  pass  into  type  2  curves.  We  show  the  ones 
for  mQ  =  6.0,  and  6.5  (they  lie  almost  on  top  of  each  other).  In  this 
tl  <  t3  <  t2  case  we  pass  directly  from  type  1  to  type  3  at  the  lower  limit  of 
the  continuum.  The  weak  saddle  point  (not  shown  in  Fig. (9))  that  lies  far  to 
the  right  of  the  node  is  not  reached  by  any  trajectory  in  this  case.  This  is 
due  to  the  behavior  of  the  F  =  0  curve  and  can  be  understood  looking  at 
Figs . ( 10)  and  ( 11 )  . 

In  Fig. (10)  the  curve  F  =  0  is  shown  for  e  =  0.0025  and  m0  =  2.5  . 
There  are  no  singular  points.  Also  shown  is  the  type  1  trajectory  for 
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T 


«0  =  2.S  .  Note  how  the  type  1  curve  starts  with  a  negative  slope  and  goes 
under  the  F  =  0  curve.  Then  it  rises  up  to  the  F  =  0  curve  and  looks  like  it 
would  become  a  type  3  curve  if  only  a  sink  were  available.  The  type  1  curve 
then  follows  the  F  =  0  curve  very  closely  until  it  finally  gets  above  it  and 
dips  down  as  a  type  1  trajectory  is  supposed  to  do.  In  Fig. (11)  we  show  the 
F  =  0  curve  for  tnQ  =  3.5,  a  little  bigger  than  in  Fig. (10).  Now  there  are  both 
a  sink  and  a  saddle  point.  What  was  a  type  1  trajectory  curve  for  mo  =  2.5  in 
Fig. (10)  has  evolved  directly  into  a  type  3  that  closely  follows  the  F  =  0 
curve  until  it  hits  the  stable  node.  This  explains  the  1 1 < 1 3  bracketing  in 
the  lower  part  of  the  continuum.  The  trajectories  are  trapped  in  the  sink  and 
cannot  approach  the  region  near  the  saddle  point.  In  this  situation  the  saddle 
point  has  little  effect  on  the  relevant  flow  regions  of  the  phase  flow. 

SUMMARY  OF  RESEARCH  RESULTS 

The  discussion  above,  illustrated  by  Figs.(l)  -  (11)  covers  all  the 
qualitative  features  that  were  encountered.  Fig. (12)  shows  the  curve  for  the 
transonic  values  of  m0  vs.  e  for  the  parameters  B  =  12.5,  G  =  3.0,  and 
RV  =  0.41  .  These  parameter  values,  as  well  as  the  others  used  here  are 
related  to  examples  given  in  Ref.l.  We  only  consider  the  detonation  region, 
which  lies  above  the  line  mQ  =  1  .  The  detonation  part  of  the  e  -  mQ  plane  is 
divided  into  3  regions  ,  labeled  1,  2,  and  3.  The  label  indicates  what  type  of 
trajectory  is  produced  by  a  pair  (e,  m0)  that  lies  in  that  region.  So,  for 
example,  if  the  pair  (e,  mQ)  lies  in  region  2  the  trajectory  associated  with 
it  is  type  2. 

The  pairs  (e,  mQ)  that  lie  on  the  curves  that  separate  the  regions 
produce  transonic  trajectories  bracketed  by  the  two  adjacent  types  of 
trajectories.  Therefore  the  curves  give  the  dependence  of  the  speed  of  steady 
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detonation  waves  vs.tbe  loss  parameter.  Region  3  is  the  continuum,  so  any 
point  in  that  region  gives  a  transonic  trajectory  of  type  3.  Fig. (13)  has  B 
and  RV  half  as  large  as  in  Fig. (12).  'Hits  represents  an  explosive  with  twice 
the  heat  of  reaction  as  the  first  one  but  the  same  activation  temperature. 
This  is  clearly  a  stronger  explosive  with  larger  detonation  speeds  but  no 
qualitative  differences.  For  both  of  these  explosives  there  is  a  lack  of  the 
weak  roots  below  the  continuum  so  region  3  is  bounded  below  by  region  1  (or, 
of  course  by  the  mQ  =  1  line). 

Fig. (14)  shows  a  case  with  the  heat  of  reaction  again  doubled  (RV  is  one 
fourth  of  the  first  example),  but  whose  activation  energy  is  twice  as  large  as 
the  other  examples.  For  small  losses,  the  detonation  speed  is  larger  (in  all 
these  cases  it  goes  to  the  ideal  Chapman-Jougue t  value)  than  the  others.  This 
case  exhibits  the  weak  solutions  below  the  continuum  so  the  detonation  part  of 
the  e  -  mo  plane  is  now  divided  into  four  regions.  The  region  of  type  2 
below  region  3  is  hardly  visible  in  Fig. (14).  Fig. (15)  is  an  expanded  view 
that  shows  the  small  region.  Another  feature  of  this  last  example  is  how 
slowly  the  upper  limit  of  region  3  goes  to  m  =  1.  In  the  previous  two 
examples  the  continuum  ends  near  e  =  0.00001  .  In  this  example  the  upper 
limit  of  the  continuum  extends  to  much  smaller  values  of  e.  The  calculation 
of  the  trajectories  in  this  neighborhood  are  the  longest  ones.  What  happens  is 
that  at  these  small  values  of  mo  the  reaction  rate  starts  very  small  and  the 
reaction  requires  a  longer  distance  to  develop.  This  is  more  pronounced  the 
larger  the  activation  temperature.  Here  wc  can  have  reaction  zone  thicknesses 
that  are  millions  of  times  longer  than  they  are  in  the  strong  detonation 
cases. 


The  physical  significance  of  model  calculations  such  as  the  one  done  in 
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this  work  is  not  clear.  In  real  explosives  the  detonation  speed  is  reduced 
when  the  losses  are  increased  in  agreement  with  the  curves  in 
Figs. (12)  -  (IS).  However.  the  limiting  point  of  detonation  failure  sets  in 
very  early  in  real  explosives,  after  the  detonation  speed  has  fallen  only 
about  10  percent.  In  these  examples,  the  limiting  point  of  detonation  failure 
occurs  after  mQ  has  dropped  by  75  percent  and  more.  Possibly  the  solutions  of 
region  3  have  an  influence  in  producing  a  much  earlier  detonation  failure.  If 
the  strong  detonation  wave  were  in  a  range  of  e  where  the  continuum  exists, 
it  might  excite  a  lot  of  the  weaker  type  3  waves  that  would  drain  away  the 
energy. 
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ABSTRACT 

Numerical  integration  of  the  steady  detonation  and  deflagration  equations 
with  a  model  expansion  loss  and  reaction  rate  can  give  four  propagation 
speeds.  The  first  two  are  a  fast  (strong)  and  a  slow  (weak)  detonation 
solutions  .  The  other  two  are  a  fast  (weak)  and  a  slow  (strong)  deflagration 
solutions.  The  two  weak  solutions  join  at  Mach  one  at  a  value  of  the  loss  that 
depends  on  the  other  parameters  of  the  problem.  These  calculations  should  be 
expanded  to  different  dependences  on  the  reaction  progress  variable. 
Approximate  analytical  results  that  confirm  the  above  numerical  calculations 
have  been  obtained.  These  results  need  to  be  pursued  because  they  could  help 
to  understand  the  influence  of  loss  terms  on  detonation  waves. 


DESCRIPTION  OF  PROPOSED  RESEARCH 

INTRODUCTION 

I  propose  to  study  the  existence  of  multiple  propagation  speeds  and 
detonation  failnre  for  steady  detonation  waves  in  the  presence  of  expansion 
losses.  The  treatment  is  based  entirely  on  the  model  developed  in  Sec.5G  of 
the  book  by  Fickett  and  Davis^.  As  discussed  in  Ref.l,  we  do  not  consider 
the  complete  two  dimensional  flow  in  a  cylindrical  explosive  charge.  This 
would  entail  the  solution  of  a  free  boundary  value  problem  of  great 
complexity.  Rather  the  equations  are  specialized  to  the  x-axis  where  the 
radial  velocity  u>  vanishes  by  symmetry.  The  divergence  du/dx  is  put  in 

as  a  function  of  x  obtained  from  some  approximation.  Much  can  be  learned  from 

the  study  of  the  remaining  equations  into  which  du/dx  is  inserted.  Their 

phase  plane  is  sufficiently  complicated  that  unexpected  behavior  is  possible. 

Ve  introduce  an  expansion  loss  parameter  defined  by 

1  dA  2  do) 

a  ■ - m - 

A  dx  u  dr 

where  u  is  the  x-component  of  fluid  velocity.  The  equations  take  the  form  of 
nozzle  flow  with  variable  cross-sectional  area  A  in  the  quasi-one-dimensional 
approximation.  The  equations  arc  valid  in  the  region  from  x  =  0,  where  the 
reaction  zone  begins,  to  the  end  of  the  reaction  zone  that  occurs  at  the  sonic 
point  u  ■  c.  As  in  Ref.l,  the  unsteady  region  connecting  the  sonic  point  to 
the  rear  boundary  condition  is  not  considered  here.  With  the  above 


restrictions,  the  equations  are 


0 


d( puA) 
dx 

do  (oR  -  ae) 
dx  "  (1  -  o2/ci) 


sod 


dX  _  R 
dx  o 

d(0  +  u2/2)/dx  =  0 


where 


H 


_ 

(r-Dp 


Xq 


with  all  the  quantities  defined  as  usual  in  Ref.  1.  These  equations  can  be  put 

into  a  for®  where  the  only  variables  are  the  Mach  number  squared, 

2  2 

m  -  u  /c  ,  and  the  reaction  progress  variable,  X.  These  equations  are 


2(l-m)dm/dx  d(^Xjq- )/dx  2e 

m(l+ym) (2+(y-l)m)  ^Xiqi+h1+u2/2  1+ym 


(1) 


and 


2 

c 


2  ( y — 1 ) 

2  ♦(y-l)m 


+  hj  +  Uj/2) 


(2) 


We  are  in  the  frame  of  the  reaction  front.  The  unreacted  explosive  fluid 
approaches  from  x  <  0  with  u  velocity  D  to  be  determined.  D  will  be  obtained 
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>■!>  .  . 


in 

terms  of 

the 

Mach  number 

squared  siq  = 

D2/ 

2 

cq  at  x  : 

-  <f. 

We  can 

discuss  both 

t&e 

detonat ion 

And 

de  f 1 agra t ion 

solutions  with 

Eqs. (1) 

and  (2). 

For 

the  detonation 

case  there 

is 

a  shock  at  x  = 

0. 

The  Mach 

number 

squared 

■l 

■  (Uj/Cj)2 

at 

x  =  0+ 

is 

related  to 

"0 

by  the 

usual 

Rugoniot 

2  * 

relations 

2  +  (y-l)m0 

1  2rm0  -  (r-D 

and 

2  2  (2r®o-(r-l))(2+(y-l)m0) 

1  0  mofy+l) 

For  the  deflagration  solutions  there  is  no  shock  at  z  =  0,  so  «•  ,  and 

c0  "  cl* 

Several  authors  have  considered  this  model  with  different  reaction  rates 

3 

and  expanson  terms  e.  Veckcn  considered  the  case  &  »  constant  and 
obtained  partial  results  on  multiple  solutions  and  detonation  failure. 
Tsuge^  et  al.,  and  Fuj iwara  and  Tsuge^  have  integrated  the  nozzle 
equations  numerically  for  a  hydrogen/oxygen  detonation  with  a  good  model  for 
the  lateral  expansion.  They  have  shown  that  sufficiently  great  expansion 
losses  lead  to  detonation  failure.  They  also  obtain  two  detonation  speeds  that 
join  at  the  failure  diameter.  For  small  losses  the  fast  root  goes  to  the 
Chapman-Jouguet  value.  They  conjectured  that  the  slow  root  went  to  Mach  one  as 
the  losses  decreased  to  zero  (charge  diameter — >  °°)  .  They  were  not  able  to 


verify  it  due  to  the  large  amounts  of  computer  time  required. 


PREVIOPS  WORK  BY  THE  PRINCIPAL  JN\T.STHWUI'OR 

I  have  obtained  preliminary  result*;  that  show  the  existence  of  2 
detonation  roots  and  2  deflagration  roots  together  with  detonation  and 
deflagration  failure.  Spec  if ica 1 1> ,  F.qs.(l)  and  (2)  have  been  integrated 
ntuaerically  with  e  =  constant  and  a  rude!  react  icr.  rate 

R  “  (c/L)  (Tt/T  ♦  l)  CXp(-Ta/T)  (1  -  X)2  ,  (3) 

where  L  it  a  reaction  length  and  T#  is  an  activation  temperature.  This  form 

for  R  was  chosen  partly  to  compare  the  numerical  results  with  the  analytical 

2 

approximation  discussed  below.  Using  T  =  M^c  /R0y  .  where  Mw  is 

molecular  weight,  and  RQ  it  the  gas  constant, we  write 


Ta 

—  -  B 


2  ♦  (y-l)m 
(y+1) ( A+h) 


where 


2  2  (2+(y-l)mn) 

h  -  (Hj+Uj/lJ/q  -  c0  7—7- 


B  -  rR0T.(r+l)/(2qMw(y-l)) 


The  phase  plane  of  Eqs.(l)  and  (2)  can  look  very  different  depending  on  the 
values  of  the  parameters  and  on  the  value  of  diq.  A  typical  case  is  shown  in 
Fig.l.  Lot  Q  be  the  right  hand  side  of  Eq.(l).  The  curves  Q  ~  0  and  m  =  1  are 
shown  in  Fig.l.  For  a  given  value  of  ni q  the  phase  orbit  starts  at  X  -  0, 
m  =  and  crosses  the  line  Q  =  0  this  is  an  orbit  of  type  1.  For  another 

value  of  biq  the  phase  orbit  reaches  the  line  m  -  1  with  a  vertical  slope 
this  la  an  orbit  of  type  2.  The  above  values  of  Oq  bracket  the  correct 
value.  The  correct  solution  is  the  one  where  dm/dx  remains  finite  at  m  -  1  so 
that  it  can  be  joined  to  the  rear  flow  behind  the  sonic  point.  The  correct 
orbit  shown  in  Fig.  1  passes  through  the  critical  point  at  m  -  1,  0-0.  The 
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problem  is  not  an  initial  value  problem  but  rather  a  boundary  value  problem 
that  yields  the  value  for  ihq.  Its  form  is  suitable  for  a  Nonlinear  Shooting 
Method^.  We  start  with  a  very  low  value  for  ,  say  0.01,  get  that  orbit, 
and  increase  Bq^  When  the  type  of  the  orbit  changes  we  have  bracketed  the 
desired  value  of  Bq  .  This  procedure  is  done  for  each  value  of  e. 

2 

With  fairly  realistic  values  for  the  parameters  Cg/q  ”  0.41, 
Y  =*  3 ,  and  B  =  12.5.  the  results  obtained  for  Bq  as  a  function  of  e  are 
shown  in  Fig.  2  where  loguiQ  is  plotted  versus  log c.  The  curve  separates 
the  regions  of  type  1  and  typo  2  in  the  logmQ-loge  plane.  Branch  A  is  the 
fast  (strong)  detonation  branch.  It  joins  the  slow  (weak)  detonation  branch  at 
Sj  ~  0.0246  .  For  e  >  e ^ ,  there  are  no  steady  solutions.  This  is 

the  value  for  detonation  failure.  The  slow  detonation  branch  ends  at  Mach  one 

at  the  value  ~  10  ^  (the  possible  singularity  at  idq  «=  1  is  not  as 

dangerous  as  it  alight  appear).  There  it  joins  the  fast  (weak)  deflagration 
branch  C.  Branch  C  joins  the  alow  (strong)  deflagration  branch  D  at 
C2  ~  2.1x10  .  This  is  the  deflagration  failure  value.  The  bend  in  the 

-7 

fast  deflagration  branch  occurs  at  ~  0.97x10  .  It  is  clear  that 
the  slow  detonation  branch  docs  not  go  to  Mach  one  at  e  =  0,  but  at 

e  ■  8j.  In  Fig.  2  one  can  see  that  for  e  <  there  are  the  two 

nornal  roots,  a  deflagration  and  a  detonation  root.  For  <  a  < 

there  are  four  roots,  while  for  (  c  (  Cj  there  are  the  two 

detonation  roots. 

The  reaction  zone  thickness  varies  widely  as  Bq  and  e  are  varied.  It 


ranges  from  the  order  of  20L  to  the  order  of  100, 00L.  For  this  reason  some 
calculations  take  Buch  more  CPU  time  than  others.  This  is  because  the  computer 
prograa  integrates  as  far  as  needed  until  it  reaches  either  □  <  0  (type  1)  or 


»  >  1  (type  2).  A  full  run  of  m0  values  for  a  particular  value  of  e  may 
take  1$  ■inotes  of  CPU  time  without  using  double  precision.  The  CPU  time  is  as 

7 

smell  es  this  thanks  to  the  use  of  a  Rungc-Kutta-Fehlberg  method  with  its 
variable  step  size  and  error  control  to  integrate  the  differential  equations. 
Fortunately  when  the  reaction  thickness  is  very  large  the  solutions  vary 
slowly  so  the  method  can  take  large  step  sizes,  say  SOL.  For  the  small 
reaction  zones,  however,  the  step  sizes  arc  much  smaller. 


PROPOSED  WORK 

I  propose  to  continue  the  numerical  calculations  just  described.  It  is 
necessary  to  use  double  precision  in  at  least  some  ranges  to  be  fairly  sure  of 
the  results. ..Also  one  should  use  a  variety  of  parameter  values  and  reaction 
rates  R.  This  would  show  how  a  different  A  dependence  in  R  affects  the 
results. 


I  would  also  like  to  pursue  obtaining  qualitative  results  that  help  to 
understand  the  numerical  ones.  For  example.  1  have  done  a  problem  with 


dX  R. 

_  „  (1_x) 

ax  u# 


where  R#  is  the  value  (constant  in  x)  of  the  R  given  in  Eq.(3)  at  the  sonic 
point.  and  u(  is  the  value  (also  constant  in  x)  of  u  at  the  sonic  point. 
With  a  further  approximation  for  the  x  dependence  of  t,  I  have  been  able  to 
integrate  Eqs.(l)  and  (2)  exactly.  The  results  are  two  nonlinear  algebraic 
equations  whose  simultaneous  solution  yields  mg.  The  results  give  a  curve 
very  much  like  that  of  Fig.  2.  Furthermore  the  algebraic  equations  obtained 
give  a  good  graphical  reason  for  why  there  should  be  4  roots. 
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PLAN  OF  WORK  AND  BUDGET  JUSTIFICATION 

I  plan  to  work  full  time  on  the  above  problems  from  July  1  to  August  15. 

1982.  and  froa  May  15  to  June  30,  1983.  This  constitutes  3  months  of  summer 

*■ 

work  during  the  one  year  period  of  this  proposal.  The  University  of  Miami  hat 
a  fast  Univac  1100/81  main  frame  computer  and  many  terminals  on  campus. 
Daring  the  day,  runs  that  request  sizable  amounts  of  CPU  time  are  given  a  low 

priority  so  the  tnrn  around  time  can  be  quite  slow.  The  runs  could  be 

submitted  in  batch  mode  but  that  destroys  the  advantages  of  demand  mode  where 
one  can  submit  a  run.  inspect  the  results,  vary  the  step  size,  and  so  on.  I 
believe  these  conditions  justify  my  request  for  a  computer  terminal  and  a 
modem  to  install  in  my  home.  This  way  I  can  submit  runs  at  night  when  the  load 
on  the  Univac  is  much  less  and  the  tnrn  around  time  is  very  short. 

STATEMENT  OF  OTHER  SUPPORT 

The  previous  work  by  the  P.I.  on  detonation  problems  was  done  mostly 
under  AFOSR  minigrant  No.AF0SR80-G135 .  That  minigrant  expired  in  the  summer 
of  1981.  The  P.I.  does  not  have  any  other  research  support  at  this  time. 
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BUDGE  i 


July  l,  1982  -  June  '  '} o i 


1.  SALARIES  AND  WAGES 

Principal  Investigator  -  Dr.  M.  Huerta 
Full  time  l|  months  summer,  102-2 
Full  time  l|  months  summer,  l')53 

Fringe  Benefits  (a  19.  3 5 

TOTAL  SALARIES  AND  WAGES 

2.  PERMANENT  EQUIPMENT 

Hazeltine  Esprit  Compute r  Terminal 
with  Koch  modifications 

1200  baud  telephone  modem 


5, 525. 
5,  Vo 7. 
11, 4?2. 

2,  224. 


13, 716. 


800. 

600. 


TOTAL  EQUIPMENT 


1, 400. 


3.  EXPENDABLE  SUPPLIES  AND  MATERIALS 


200. 


4.  TRAVEL  (DOMESTIC) 

One  scientific  meeting  and  one  visit 
to  another  institution; 

Estimated  r.t.  air  fare  fe-  $350/trip 

2  x  $350 

3  days  per  diem/trip  ('•  $75/day 

2  x  3  x  $75 

Local  transportation  (“  $25 /trip 
2  x  $25 

TOTAL  TRAVLL 

5.  PUBLICATION  COSTS 


700. 

450. 

50. 


1, 200 

700 
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July  ),  1  *M.’  -  June  vi',  1  ' 1  } 


( <  ontinm.  !} 


6.  OTHER  COSTS 

— 

Computer  Time 


1,  bOO. 


Long  Distance  Telephone  _ 50, 

TOTAL  OTHER  COSTS  1,  550. 

7.  TOTAL  DIRECT  COSTS  18,  766. 

8.  INDIRECT  COSTS  (see  note  below) 

60%  Modified  Total  Direct  Costs  (.  60  x  $17,  366)  _ 10,  420. 

9.  TOTAL  AMOUNT  REQUESTED  29,  186. 


NOTE:  The  Coral  Gables  campus  negotiated  indirect  cost  rate  for  the 

University  of  Miami  for  fiscal  year  6/1/81  to  5/31/ 82  is  94.  0% 
of  Modified  Total  Direct  Costs  (MTDC).  The  University  is 
willing  to  accept  an  indirect  cost  rate  that  is  lower  than  the 
negotiated  rate.  However,  if  the  proposal  is  awarded,  the 
University  retains  the  right  to  exercise  any  cost  clause  or 
cost  regulation  made  a  part  of  the  signed  agreement. 
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Name:  Manuel  A.  Huerta  Tl  tic :  Associate  Professor  of  Physi 


bate:  November  30,  1980 


Professional  Exoerience: 


Associate  Professor  of  Physics 
University  of  Miami,  Coral  Gables,  Florida 
September  1978 

Program  Associate  (on  leave  from  University  of  Miami) 

U.S.  -  Latin  America  Cooperative  Science  Program 
National  Science  Foundation 
May  1977  -  August  1973 

Associate  Professor  of  Physics 
University  of  Miami,  Coral  Gables,  Flori  la 
May  1975  -  (Tenured  in  May  1976) 

Assistant  Professor  of  Physics 
University  of  Miami,’  Coral  Cables 
September  1972  -  May  1975 

Associate  F.esearch  Scientist,  Magneto  Fluid  Dynamics  Division 
and  Adjunct  Professor  of  Mathematics 

Courant  Institute  of  Mathematical  Sciences,  New  York  University 
September  1970  -  August  1972 

Instructor  of  Physics 
University  of  Miami,  Coral  Gables 
September  1966  -  August  1970 


Education : 

Ph.D.  Physics,  January  1970:  University  of  Miami,  Coral  Gables,  Florid 
Grade  Point  Average:  4.0 

Dissertation:  Approach  to  Fguilibriun  of  Infinite  Chains  of  Coupled 

Harmonic  Oscillators 


M.S.  Physics,  January  1967;  University  of  Miami,  Coral  Gables,  Florida 
Grade  Point  Average:  4.0 
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B.S.  Electrical  Engineer ing,  June  19G5:  California  Institute  of 
Technology,  Pasadena,  California 
Grade  Point  Average:  3.3 
Graduated  with  Honors 

Enrolled  at  University  of  Miami  in  Electrical  Engineering  from 
September  1961  to  August  1963  then  transferred  to  the  California 
Institute  of  Technology. 

Grade  Point  Averages  3.8 


Scholarships : 

NASA  Fellowship  in  Physics,  University  of  Miami,  Coral  Gables 
September  1965  -  August  1966. 

Graduate  Assistant  in  Physics,  Center  for  Theoretical  Studies, 
University  of  Miami,  Coral  Gables 
Summer  of  1965. 

Undergraduate  Scholarship,  California  Institute  of  Technology 
1963  -  1965. 

Undergraduate  Scholarship,  University  of  Miami,  Coral  Gables 

1962  -  1963. 


Honors  as  Student: 


California  Institute  of  Technology's  David  Joseph  McPherson  Prize. 
Awarded  upon  graduation  to  the  graduating  senior  in  engineering 
who  best  exemplified  excellence  in  scholarship. 


Membership  in  Societies.- 

Tau  Beta  Pi,  Phi  Kappa  Phi,  Sigma  Xi,  American  Physical  Society, 
American  Mathematical  Society. 


Research  Fields  of  Interest: 


Plasma  Physics,  Statistical  Mechanics,  Solid  State  Physics, 
Acoustic  Propagation,  Nonlinear  Partial  Differential  Equations. 


Professional  Consulting: 

September  1973  -  August  1974:  Consultant  to  the  Atomic  Energy 
Commission  through  the  Magneto  Fluid  Dynamics  Division  of  the 
Courant  Institute  of  Mathematical  Sciences  of  New  York  University. 
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Participation  in  Research  Cent. rn cts : 

September  1973  -  August  197-1:  Principal  Investigator  in  project 
entitled  "Research  in  Acoustic  Propagation  and  Reflection  from 
Internal  Wavts"  done  under  a  MO  A  A  Grant  Number  04--4-022-9. 

September  1974  -  August  1975:  Renewal  of  above  grant. 

September  1975  -  August  197G:  Renewal  of  above  grant. 

August  1979  -  July  1980:  Principal  Investigator  in  Contract 
Number  AFOSR  78-3663  entitled  "The  Stability  and  Dynamics  of 
Elastic  Structures  and  Fluid  Flows." 


July  1980  -  July  1981:  Principal  Investigator  in  Grant  AFOSR  80-0135 
entitled,  "Analytical  Studies  of  Nonidcal  Explosives" 


Juried  Journal  Articles 
of 

Manuel  A.  Huerta 


1.  With  T.N»  Veziroylu,  Thermal  Conductance  of  Two  Dimensional 
Eccentric  Constrictions;  NASA  STAR  6,  20,  October  23,  1963; 
STAR  Index  N-68-32974'. 

2.  With  H.S.  Robertson,  Entropy,  Information  Theory  and  the 

Approach  to  Equilibrium  of  Coupled  Harmonic  Oscillator 
Systems’;  J~.  Stat~.  Phys .  393  (1969)  . 

3.  With  H.S.  Robertson,  Approach  to  Equilibrium  of  Coupled 
Harmonically  Bound  Systems;  Phys.  Rev.  Letters  23,  825 
(19  6  971 

4.  With  H.S.  Robertson,  Information  Theory  and  the  Approach 
to  Equilibrium;  Am.  J.  Phys.  3_8  619  (1970). 

5.  With  H.S.  Robertson,  Entropy  Oscillations  and  the  H  Theorem 
for  Finite  Segments  of  Infinite  Coupled  Harmonic  Oscillator 
Chains;  Pure  and  Appl.  Chem.  22_,  4~09  (1970)  . 

6.  With  H.S.  Robertson,  Approach  to  Equilibrium  of  Coupled 
Harmonic  Oscillator  Systems  II;  J.  Stat.  Phys.  3,  171 
(1971)1 

7.  With  H.S.  Robertson  and  J.C.  Nearing,  Exact  Equilibration 
of  Harmonically  Bound  Osfcillator  Chains;  J.  Math.  Phys. 

12,  2305  (1971).  — 

8.  With  G.  Gonzalez,  Field  Expressions  for  a  Circular  Loop 
Antenna  in  Terms  of  a  New  Set  of  Functions;  J.  Appl.  Phys. 

43,  3975  (197277 

9.  With  G.  Gonzalez,  Fresnel  Region  Fields  Produces  by  Circular 
Loop  and  Annular  Slot  Antennas;  J.  Appl.  Phys .  44,3500  ( 1973) 

10.  With  T.N.  Veziroglu  and  S.  Kakac,  Exact  Solutions  for  Thermal 
Conductances  of  Two  and  Three  Dimensional  Contacts;  Thermal 
Conductivity  14  (1976),  Edited  by  P.G.  K1 omens  and  T.K.  Chu, 
Plenum  Publishing  Corp.,  Now  York,  New  York. 

11.  With  G.C.  Alcxandrakis  and  O.L.S.  Lieu,  Theory  of  Nonlinear 
Phenomena  in  Ferromagnetic  Transmission  Resonance;  Phys.  Rev. 
B16,  476  (1977) . 
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12. 


13. 


With  R.A.  Goldstein  and  J.C.  Nearing.  ||^ti£^a 


TrTL  kraan  Plasma  as  a  Nonlinear  Bifurcation, 
Phys.  Fluids  *  2  2  2  31,  ( 1 9  7  9 )  . 


ry  Striations 
Phenomenon ; 


With  G.  Gonzalez,  Surface  Potentia ls_  of. 
Volume  Conductors  Excited  fay  a  Current^ 
Electronics  iT~,  213  (1979)  . 


Stratified  Spheroidal 
Dipole  Source;  int.  3. 


14.  With  J.  Magnan,  Spatial  Structures  in  Plasmas  with  Metastable 
States  as  Bifurcation  Phenomena;  Accepted  for  publication 
Phys .  Rev.  A 


15.  With  G.  Gonzalez,  The  Surface  Potentials  Produced  by  Electric 
Sources  in  Stratified  Spherical  and  Prolate  Spheroidal  Volume 
Conductors ;  Submitted  to  IEEE  Trans.  Biom.  Eng. 
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Other  I  ubi i 1  Li  :;nj 


by 

ndiniol  A.  Hu  or 'a 


1.  With  J.C.  Nearing,  W.  15.  Pardo  and  fi.  Reused,  The _ Interactions 

of  Electromagnetic  Radiation  with  a  Forest  Environment;  sub¬ 
con  tract**lIrlder—Lrnco  1  n  LahoTatory  of  m.-iTt~  purchase  f order 
number  C-942  under  prime  contract  number  AF  19  ( 628 ) - 5167 . 
Available  as  a  Department  of  physics,  University  of  Miami 
Report,  MIAPII-70 . 10  . 

2.  With  J.C.  Nearing,  W.B.  Pardo  and  M.  Reusch,  Some  Physics  of 
a  Florida  Forest;  Department  of  Physics,  University  of  Miami 
Report  MIAPH-70 . 11 . 

3 .  Approach  to  Equilibrium  of  In finite  Chains  of  Coupled  Harmonic 
Oscillators ;  doctoral  dissertation,  Department  of  Physics, 
University  of  Miami  Report  MIAPH-69.24. 
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